TUTORIAL 1 : NVT Lennard-Jones fluid¶
| Authors: | Andrey V. Brukhno (andrey.brukhno{at}stfc.ac.uk), James Grant (r.j.grant{at}bath.ac.uk), and John Purton (john.purton{at}stfc.ac.uk) |
|---|
MC method in outline¶
In this tutorial we will perform the most straightforward Monte Carlo (MC) simulation - a simulation of a Lennard-Jones fluid in the canonical (NVT) ensemble. The triplet NVT signifies that number of particles, volume (hence, density) and temperature are kept constant during the entire length of this type of simulation.
In this case only the particles are allowed to move, and a random walk is constructed following the approach of Metropolis et al. [1]. The Metropolis algorithm iterates the following three steps that together constitute a single MC move (or attempt):
- Select a particle at random and calculate its energy \(U(\mathbf{r}_1)\).
- Displace the particle randomly to a trial position and calculate its new energy \(U(\mathbf{r}_2)\).
- Accept or reject the particle displacement according to the Metropolis acceptance rule:
How do we achieve this?¶
Ensure detailed balance: it is crucial to ensure a truly random (uniform) selection of trial particle positions! Otherwise an unwanted, hidden bias can be introduced into the MC sampling procedure, which would invalidate the obtained statistics.
The common way to move (translate) a particle is to generate a random displacement vector \(\Delta\mathbf{r}\) bound within a cube with a given side legnth \(\Delta_{max}\) (the so-called displacement parameter) and add it to the current position of the particle under trial.
Then, pick up a random number between 0 and 1 and accept or reject the move as illustrated below:
DL_MONTE Input Files¶
The purpose of this tutorial is to introduce the basic file formats that are used by DL_MONTE and then run an example MC simulation. DL_MONTE requires at least 3 standard input files (similar to DL_POLY):
FIELD – units, atoms, molecules, topology and force-field specs (fixed format)
CONFIG – initial configuration for the system, including cell vectors (fixed format)
CONTROL – simulation parameters: external conditions, options and flags (free format)
NOTE: The file names are fixed, i.e. only these file names are recognised by DL_MONTE. In some more intricate cases other input files might be needed (e.g. for restarting a simulation for continuation).
We will introduce these files and discuss each of them in turn. Full details of their structure, formats and default values can be found in the DL_MONTE manual.
To check the input files for the current tutorial, navigate to the corresponding directory:
[DL_MONTE-2_tutorial]$ cd exercises/tutorial_1
[DL_MONTE-2_tutorial]$ ls -l
NOTE: The first line in all the three files is an arbitrary title which can contain up to 80 characters (any symbol over that limit will be ignored by DL_MONTE). The titles can be different in each file.
FIELD¶
In DL_MONTE the system is abstracted into atoms and molecules: atoms are treated as point-like particles, and molecules are collections of atoms which can be moved collectively. This, along with a versatile selection of potential forms which can be combined into different force fields, allows for simulation of a wide range of systems comprised of possible combinations of free unconnected atoms (so-called atomic/ionic fields), and molecules possessing structure: both rigid and flexible.
The FIELD file defines all the atomic and molecular types present in the system, but also possibly the molecular structure (so-called topology: bonds, angles, etc). Finally, it also specifies the interatomic and (optional) external potentials. A detailed description of all the available force-field parameters and their interplay can be found in the DL_MONTE manual.
The file structure is similar to that used by DL_POLY, but there are some important differences owing to the requirements for Grand-Canonical simulations (considered later).
The FIELD file for our first tutorial is shown below
1 2 3 4 5 6 7 8 9 10 11 12 13 | Lennard-Jones, 2.5*sigma cut-off, sigma = 1 Angstrom, epsilon = 1 internal unit = 10 J/N_A
CUTOFF 2.5 # cut-off radius (must be the 2-nd line - after the title)
UNITS internal # energy unit (must be the 3-rd line - after CUTOFF)
NCONFIGS 1 # number of configurations in CONFIG file (must be the 4-th line)
ATOM TYPES 1 # number of atom types or elements (must be the 5-th line)
LJ core 1.0 0.0 # mass and charge of the atom(s)
MOLECULE TYPES 1 # number of molecular types
lj # molecule name (case sensitive)
MAXATOM 512 # max number of atoms within current molecular type
FINISH # closing the molecular type specs
VDW 1 # number of short-ranged Van der Waals potentials
LJ core LJ core lj 1.0 1.0 # atom names and types for a VDW pair, LJ epsilon and sigma
CLOSE # closing FIELD
|
The important points to note are:
the cutoff = 2.5 for the short-ranged VDW interactions is in Angstroms
the energy unit used in this case is internal = 10 J/mol [or 0.01 kJ/mol],
the CONFIG file contains NCONFIGS = 1 configuration.
the only ATOM type present has the name LJ and the type core, having mass = 1.0 and charge = 0.0.
All the atoms in the system are contained within a single molecule called lj. The number of atoms in this molecule is restricted to 512 at most. Note that the CONFIG file cannot contain more than that number of atoms within that type of a molecule. As a matter of fact, for an NVT simulation where the number of particles cannot change, the CONFIG file must contain the total number of atoms that is specified in the FIELD file!
The atoms interact through a Lennard-Jones interaction potential, which is specified by the keyword lj, and the LJ parameters are set to: epsilon = 1.0 and sigma = 1.0.
CONFIG¶
The CONFIG file contains the starting configuration specification. It defines the cell type and geometry, sets the cell vectors (in the matrix form), the actual and maximum numbers of molecules, atoms within each molecule, and also the coordinates of all the atoms.
Below we illustrate the structure of the CONFIG file with a sample of the configuration provided for this tutorial.
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 | Exercise 1: LJ fluid in NVT ensemble
0 0
11.7452 0.00000 0.00000
0.00000 11.7452 0.00000
0.00000 0.00000 11.7452
NUMMOL 1 1
MOLECULE lj 512 512
LJ core
0 0 0
LJ core
0 0 0.125
LJ core
0 0 0.25
LJ core
0 0 0.375
|
In line 2 the first number is an integer switch operating in the same way as levcfg in DL_POLY. It should be set to zero if the CONFIG file only contains the particle coordinates (which is sufficient for starting an MC simulation). However, if non-zero, this flag allows DL_MONTE to read in CONFIG files originating from MD (DL_POLY) simulations and containing (x,y,z) components of the particle velocities (levcfg = 1) or velocities and forces (levcfg = 2) which would appear on the lines directly following the particle coordinates. Note, however, that DL_POLY CONFIG files cannot be used directly, without certain modifications, as input for DL_MONTE simulations.
The second integer flag determines the convention used for the particle coordinates: 0 for fractional representaion, and non-zero for Cartesian: 1 == cubic, 2 == orthorhombic, 3 == palleliped. These are followed by the cell matrix.
In our case the CONFIG file contains fractional coordinates and, since the cell matrix is diagonal and all its diagonal elements are equal, the simulation cell is cubic.
The line starting with NUMMOL keyword specifies the total number of molecules present, followed by the maximum number of molecules allowed for each molecular type defined in the FIELD file. In this case we have only one molecular type lj, and there is only one molecule of this type in the CONFIG file. The molecule comprises 512 atoms (LJ particles) which equals the maximum number.
Finally, the atoms’ specs are read in. The atom data follows and for each particle, its name (usually chemical symbol) and type (core, semi and metal that can be abbreviate to c, s, and m) are specified. On the next line these are followed by the coordinates in a format consistent with that given on line 2.
CONTROL¶
The CONTROL provides directives to DL_MONTE how to undertake the calculations and switches on or off functionality. The CONTROL file in this example is:
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 | NVT simulation of Lennard-Jones fluid
#use ortho
finish
seeds 12 34 56 78 # Seed RNG seeds explicitly to the default
nbrlist auto # Use a neighbour list to speed up energy calculations
maxnonbondnbrs 512 # Maximum number of neighbours in neighbour list
temperature 1.4283461511745 # T* = 1.1876 = T(K) * BOLTZMAN (see constants_module.f90)
steps 10000 # Number of moves to perform in simulation
equilibration 0 # Equilibration period: statistics are gathered after this period
print 1000 # Print statistics every 'print' moves
stack 1000 # Size of blocks for block averaging to obtain statistics
sample coord 10000 # How often to print configurations to ARCHIVE.000
revconformat dlmonte # REVCON file is in DL_MONTE CONFIG format
archiveformat dlpoly4 # ARCHIVE.000/HISTORY.000/TRAJECTORY.000 format
# In this case: HISTORY.000 in DLPOLY4 style
move atom 1 100 # Move atoms with a weight of 100
LJ core
start
|
The first line is the title and the second contains the keyword finish. We will see later that there a number of directives in DL_MONTE that can only occur before this directive which must be present in the CONTROL file.
Seeds followed by a series of 4 integers provides a reproducible seed, otheriwse ranseed generates a random seed from the system clock at initialisation.
The diretives nbrlist and maxnonbondnbrs control the size and administration of the neighbourlist used by DL_MONTE to optimise performance, and are explained further in one of the exercises.
Temperature (K) is self explanatory while steps is the length of the simulation in attempted moves. NOTE: Providing that LJ epsilon = 1.0 in FIELD, the reduced LJ temperature T* = T(K)*R, where R = BOLTZMAN (see constants_module.f90 or beta = 1/RT(K) in OUTPUT.000). In general, T* = RT(K)/epsilon = 1/(epsilon*beta), where epsilon is in internal energy units.
Equilibration etc control the detail of how data is output from DL_MONTE.
The revconformat dlmonte instruction describes the format of the output file REVCON.000 which contains the final configuration of the simulation and can be used for continuing a simulation, dlpoly2, dlpoly4 are other options.
archiveformat dlpoly4 describes the format of the trajectory file here it will be HISTORY.000, equivalent of the HISTORY in DLPOLY4. The dlpoly{2,4} formats are readable by common visualisation packages such as vmd. Full options are detailed in the manual.
The directive move atoms is where we begin to touch on the fundamental control of the simulation. The key feature here is that DL_MONTE will not do anything unless told to do so (N.B. While this gives DL_MONTE great flexibility it means also means that it may be possible to ask DL_MONTE to perform ill-defined calculations). In this simple NVT ensemble only the particles move, thus move atoms tells DL_MONTE to move 1 atom type, with a weight of 100. The line(s) following this detail the atom and type.
Finally the start directive ends the CONTROL file and instructs DL_MONTE to start the simulation.
We are now ready to run DL_MONTE:
[tutorial_1]$ DL_MONTE.SRL.X &
When running DL_MONTE locally (i.e. in the “interactive” mode), to check if the simulation is still running, we can use the top command:
[tutorial_1]$ top
which produces a list of the active processes. If DL_MONTE.SRL.X is still running its name will be found in the list.
Press q key to quit the interactive top output.
Examining the output files¶
A successful DL_MONTE calculation will produce a number of output files:
* OUTPUT.000 - details of the simulation, statistics, running time, or errors if the calculation failed.
* REVCON.000 - the final configuration in the format specified
* PTFILE.000 - statistics though will eventually be deprecated in favour of
* YAMLDAT.000 - statistics in the yaml format
* ARCHIVE.000/HISTORY.000/TRAJECTORY.000 - the trajectory in the specified format
For analysis we will typically process the YAMLDAT.000 and visualise the trajectory files. However for understanding how the simulation proceeds it is useful to have some familiarity with the OUTPUT file. The file begins with a header detailing the version, authors and suggested citations, followed by the brief summary of details of the simulation as specified in the input files. A section headed simulation parameters then specifies all parameter values that will be used within the simulation.
The final step before starting the calculation is to determine the initial energy of the system and the details of this are printed in a block:
--------------------------------------------------
initial energies
--------------------------------------------------
break down of energies for box: 1
total energy -0.7037212307E+03
reciprocal space coulomb 0.0000000000E+00
real space coulomb 0.0000000000E+00
external mfa coulomb 0.0000000000E+00
nonbonded two body (vdw) -0.7037212307E+03
bonded two body (pair) 0.0000000000E+00
nonbonded three body 0.0000000000E+00
bonded three body (angle) 0.0000000000E+00
bonded four body (angle) 0.0000000000E+00
many body energy 0.0000000000E+00
external potential energy 0.0000000000E+00
total virial 0.0000000000E+00
volume 0.1620247087E+04
This is followed by a partial breakdown per molecule type in the system and the time taken to initialise the calculation. There after every print steps as specified in the CONTROL file a block is printed to the OUTPUT.000 file:
Iteration 10000 - elapsed time: 0 h : 0 m : 0.301 s
====================================================================================================
step en-total h-total coul-rcp coul-real
step en-vdw en-three en-pair en-angle
step en-four en-many en-external en-extMFA
step volume cell-a cell-b cell-c
step alpha beta gamma
r-av en-total h-total coul-rcp coul-real
r-av en-vdw en-three en-pair en-angle
r-av en-four en-many en-external en-extMFA
r-av volume cell-a cell-b cell-c
r-av alpha beta gamma
----------------------------------------------------------------------------------------------------
10000 -0.1006733261E+04 -0.3030120301E+03 0.0000000000E+00 0.0000000000E+00
-0.1006733261E+04 0.0000000000E+00 0.0000000000E+00 0.0000000000E+00
0.0000000000E+00 0.0000000000E+00 0.0000000000E+00 0.0000000000E+00
0.1620247087E+04 0.1174520000E+02 0.1174520000E+02 0.1174520000E+02
0.9000000000E+02 0.9000000000E+02 0.9000000000E+02
-0.1014992775E+04 -0.3112715438E+03 0.0000000000E+00 0.0000000000E+00
-0.1014992775E+04 0.0000000000E+00 0.0000000000E+00 0.0000000000E+00
0.0000000000E+00 0.0000000000E+00 0.0000000000E+00 0.0000000000E+00
0.1620247087E+04 0.1174520000E+02 0.1174520000E+02 0.1174520000E+02
0.9000000000E+02 0.9000000000E+02 0.9000000000E+02
LJ c 512.0000 512.0000
lj 1.0000 1.0000
----------------------------------------------------------------------------------------------------
equilibration period ended at step 10000
----------------------------------------------------------------------------------------------------
This specifies the same breakdown in a tabulated form, including the iteration number and time taken. By using the command:
[tutorial_1]$ grep time OUTPUT.000
you can quickly gauge the progress of your calculation.
DL_MONTE runs a check of the system at regular intervals to verify that the system is behaving correctly. This typically involves calculating the energy from scratch and comaring it with the running total, the result is printed to the OUTPUT.000 with a line looking like:
[tutorial_1]$ grep U_recalc OUTPUT.000
Workgroup 0, box 1 check: U_recalc - U_accum = 0.35129E-10 0.29580E-10 0.57773E-13 -0.34894E-13 (internal, kT, kT/atom, dU/U)
The final value is the relative energy difference which should be of the order of the working precision, typically ending E-13 or E-14.
Finally at the end of the simulation a summary block detailing the average values during the simulation (excluding the first equilbration steps) and their fluctations, the final energies, and ‘processing data’ detailing time, move data and final parameters:
----------------------------------------------------------------------------
processing data
----------------------------------------------------------------------------
total no of attempted & empty atom moves : 110000 0
successful no of atom moves & acceptance : 42987 0.39079091
displacement (Angstroms) for LJ c : 0.4160
pure simulation time (excl. init.): 0 h : 0 m : 3.492 s
total elapsed time (incl. init.): 0 h : 0 m : 3.498 s
normal exit
Upon a successful completion of a simulation run, the OUTPUT.000 file will end with the string normal exit. To quickly check for this try:
[tutorial_1]$ grep "normal exit" OUTPUT.000
Basic analysis of the output¶
To examine the energy evolution along the MC time, first run the script:
[tutorial_1]$ strip_yaml.sh energy
and then plot the energy vs MC time:
[tutorial_1]$ gnuplot
gnuplot> plot './energy.dat' u 1:2 w l t "Energy(MC step)"
To visualise the trajectory we can use vmd:
[tutorial_1]$ vmd
To visualise the trajectory:
File
New molecule
Determine file type
-> Browse and Select HISTORY.000
Select the correct format:
-> DLPOLY V3 History
Note that in a newer version of VMD it is called DLPOLY-4
Questions to ask yourself:¶
- What was the initial configuration, and what did happen during simulation?
- Did the system reach equilibrium?
- Can we improve on the stats? How?
Extra exercises to try and do¶
Now try the extension exercises to learn more about funcitonality within DL_MONTE and to optimise your calculation:
With each exercise we recommend you create a copy of the inputs in a sub-directory:
[tutorial_1]$ mkdir ex1
[tutorial_1]$ cp CONFIG CONTROL FIELD ex1
[tutorial_1]$ cd ex1
Or move on to TUTORIAL 2 : NPT Lennard-Jones fluid and NPT simulations.
Footnotes
| [1] |
|